Chapter 4 Community composition
4.1 Taxonomy overview
4.1.1 Stacked barplot
# Merge data frames based on sample
merged_data<-genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
filter(count > 0) #filter 0 counts
# Create an interaction variable for time_point and sample
merged_data$interaction_var <- interaction(merged_data$sample, merged_data$time_point)
ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
facet_nested(. ~ type, scales="free") + #facet per day and treatment
guides(fill = guide_legend(ncol = 1)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.title.x = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
labs(fill="Phylum",y = "Relative abundance",x="Samples")+
scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
facet_wrap(~type, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show time_point label
4.1.2 Phylum relative abundances
phylum_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
left_join(genome_metadata, by = join_by(genome == genome)) %>%
group_by(sample,phylum) %>%
summarise(relabun=sum(count))
phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| phylum | mean | sd |
|---|---|---|
| p__Bacteroidota | 0.380806816 | 0.202118047 |
| p__Bacillota_A | 0.253813952 | 0.163904042 |
| p__Bacillota | 0.118079729 | 0.146679823 |
| p__Pseudomonadota | 0.096901291 | 0.160581847 |
| p__Campylobacterota | 0.053050348 | 0.092998130 |
| p__Verrucomicrobiota | 0.029867960 | 0.073128363 |
| p__Desulfobacterota | 0.023580475 | 0.036482484 |
| p__Chlamydiota | 0.010572704 | 0.059690269 |
| p__Fusobacteriota | 0.010482132 | 0.028311738 |
| p__Cyanobacteriota | 0.009206465 | 0.016492133 |
| p__Bacillota_C | 0.004713164 | 0.006645703 |
| p__Spirochaetota | 0.004009680 | 0.012308028 |
| p__Bacillota_B | 0.002465465 | 0.004927667 |
| p__Actinomycetota | 0.001235741 | 0.006346852 |
| p__Elusimicrobiota | 0.001214079 | 0.006501752 |
phylum_arrange <- phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun)) %>%
arrange(-mean) %>%
select(phylum) %>%
pull()
phylum_summary %>%
filter(phylum %in% phylum_arrange) %>%
mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
geom_jitter(alpha=0.5) +
theme_minimal() +
theme(legend.position="none") +
labs(y="Phylum",x="Relative abundance")
4.2 Taxonomy boxplot
4.2.1 Family
family_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,family) %>%
summarise(relabun=sum(count))
family_summary %>%
group_by(family) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| family | mean | sd |
|---|---|---|
| f__Bacteroidaceae | 2.215146e-01 | 0.1392048812 |
| f__Lachnospiraceae | 1.406192e-01 | 0.1085410432 |
| f__Tannerellaceae | 1.028077e-01 | 0.0799387840 |
| f__Helicobacteraceae | 5.258990e-02 | 0.0925721039 |
| f__Mycoplasmoidaceae | 3.694491e-02 | 0.0756423660 |
| f__Erysipelotrichaceae | 3.502391e-02 | 0.0451216980 |
| f__UBA3700 | 3.418595e-02 | 0.0557485850 |
| f__Marinifilaceae | 2.774048e-02 | 0.0271999075 |
| f__ | 2.772765e-02 | 0.0838376915 |
| f__DTU072 | 2.700317e-02 | 0.0529029477 |
| f__Rikenellaceae | 2.696900e-02 | 0.0464477971 |
| f__Enterobacteriaceae | 2.691677e-02 | 0.0918386902 |
| f__Coprobacillaceae | 2.551773e-02 | 0.0892179821 |
| f__Desulfovibrionaceae | 2.358047e-02 | 0.0364824837 |
| f__Ruminococcaceae | 1.859938e-02 | 0.0429483670 |
| f__LL51 | 1.759437e-02 | 0.0687507595 |
| f__Rhizobiaceae | 1.596807e-02 | 0.0770633474 |
| f__UBA3830 | 1.578817e-02 | 0.0454016753 |
| f__Akkermansiaceae | 1.227359e-02 | 0.0316838689 |
| f__Chlamydiaceae | 1.057270e-02 | 0.0596902690 |
| f__Fusobacteriaceae | 1.048213e-02 | 0.0283117384 |
| f__CAG-239 | 9.141331e-03 | 0.0150992357 |
| f__Enterococcaceae | 8.420523e-03 | 0.0466580889 |
| f__Gastranaerophilaceae | 7.728770e-03 | 0.0144404459 |
| f__Oscillospiraceae | 6.643218e-03 | 0.0078105983 |
| f__UBA1997 | 6.378408e-03 | 0.0309668631 |
| f__Streptococcaceae | 6.366441e-03 | 0.0342174908 |
| f__UBA1242 | 4.673359e-03 | 0.0153730061 |
| f__Brevinemataceae | 4.009680e-03 | 0.0123080282 |
| f__Acutalibacteraceae | 3.374550e-03 | 0.0109560112 |
| f__UBA660 | 3.196511e-03 | 0.0139558713 |
| f__Clostridiaceae | 2.842205e-03 | 0.0171339212 |
| f__RUG11792 | 2.817729e-03 | 0.0251127757 |
| f__Peptococcaceae | 2.465465e-03 | 0.0049276669 |
| f__MGBC116941 | 2.160169e-03 | 0.0093827415 |
| f__Acidaminococcaceae | 1.943804e-03 | 0.0050321931 |
| f__CAG-508 | 1.807978e-03 | 0.0064175379 |
| f__Anaerovoracaceae | 1.569531e-03 | 0.0036471440 |
| f__Moraxellaceae | 1.485415e-03 | 0.0097446552 |
| f__RUG14156 | 1.477695e-03 | 0.0045847831 |
| f__Staphylococcaceae | 1.361709e-03 | 0.0050863070 |
| f__Elusimicrobiaceae | 1.214079e-03 | 0.0065017516 |
| f__CAG-288 | 9.490865e-04 | 0.0060198775 |
| f__Anaerotignaceae | 8.989745e-04 | 0.0040469504 |
| f__CALVMC01 | 7.516846e-04 | 0.0043609744 |
| f__Eggerthellaceae | 6.407882e-04 | 0.0021266467 |
| f__Massilibacillaceae | 6.235073e-04 | 0.0016350480 |
| f__Mycobacteriaceae | 5.949530e-04 | 0.0060400054 |
| f__UBA1820 | 4.736030e-04 | 0.0013057353 |
| f__Arcobacteraceae | 4.604457e-04 | 0.0050324221 |
| f__CAG-274 | 4.519746e-04 | 0.0022028477 |
| f__Burkholderiaceae_C | 3.699430e-04 | 0.0048092594 |
| f__Muribaculaceae | 3.617894e-04 | 0.0009776409 |
| f__UBA932 | 3.419549e-04 | 0.0011583178 |
| f__Hepatoplasmataceae | 2.989107e-04 | 0.0038858387 |
| f__Rhodobacteraceae | 2.959092e-04 | 0.0038468195 |
| f__Weeksellaceae | 2.771627e-04 | 0.0031476408 |
| f__Eubacteriaceae | 1.646822e-04 | 0.0006729067 |
| f__Sphingobacteriaceae | 1.505775e-04 | 0.0012460017 |
| f__Devosiaceae | 1.489995e-04 | 0.0015094318 |
| f__Pumilibacteraceae | 1.277418e-04 | 0.0007646754 |
| f__WRAU01 | 9.603359e-05 | 0.0012484367 |
| f__Peptostreptococcaceae | 2.287338e-05 | 0.0002973539 |
family_arrange <- family_summary %>%
group_by(family) %>%
summarise(mean=sum(relabun)) %>%
arrange(-mean) %>%
select(family) %>%
pull()
family_summary %>%
left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
filter(family %in% family_arrange[1:20]) %>%
mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
scale_color_manual(values=phylum_colors[-8]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")
4.2.2 Genus
genus_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,genus) %>%
summarise(relabun=sum(count)) %>%
filter(genus != "g__")
genus_summary %>%
group_by(genus) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()| genus | mean | sd |
|---|---|---|
| g__Bacteroides | 1.352836e-01 | 0.0923254347 |
| g__Parabacteroides | 9.681305e-02 | 0.0801790006 |
| g__Phocaeicola | 6.935126e-02 | 0.0794715561 |
| g__Mycoplasmoides | 3.051289e-02 | 0.0753256189 |
| g__Helicobacter_J | 3.009004e-02 | 0.0594808174 |
| g__Odoribacter | 2.552307e-02 | 0.0267570773 |
| g__Roseburia | 2.384718e-02 | 0.0559439989 |
| g__NHYM01 | 2.249986e-02 | 0.0796951054 |
| g__Alistipes | 2.208351e-02 | 0.0283591693 |
| g__Coprobacillus | 2.013859e-02 | 0.0878851683 |
| g__Agrobacterium | 1.596807e-02 | 0.0770633474 |
| g__Akkermansia | 1.227359e-02 | 0.0316838689 |
| g__Fusobacterium_A | 1.038903e-02 | 0.0283170720 |
| g__Kineothrix | 8.801109e-03 | 0.0409037874 |
| g__Proteus | 8.657984e-03 | 0.0681831315 |
| g__Dielma | 8.578108e-03 | 0.0090836958 |
| g__CAG-95 | 8.109350e-03 | 0.0205019603 |
| g__JAAYNV01 | 7.994602e-03 | 0.0195404994 |
| g__UBA866 | 7.260325e-03 | 0.0293499698 |
| g__Desulfovibrio | 7.018155e-03 | 0.0211479806 |
| g__Enterococcus | 7.001300e-03 | 0.0456600395 |
| g__Ureaplasma | 6.432013e-03 | 0.0138012183 |
| g__Lactococcus | 6.366441e-03 | 0.0342174908 |
| g__Lacrimispora | 6.064761e-03 | 0.0098190433 |
| g__Parabacteroides_B | 5.994680e-03 | 0.0100174097 |
| g__CALXRO01 | 5.765728e-03 | 0.0308524398 |
| g__Citrobacter | 5.717035e-03 | 0.0334547978 |
| g__NSJ-61 | 5.560111e-03 | 0.0198839325 |
| g__Breznakia | 5.504528e-03 | 0.0237016548 |
| g__Clostridium_AQ | 5.326190e-03 | 0.0121695008 |
| g__Salmonella | 5.133583e-03 | 0.0184556360 |
| g__Bilophila | 4.983840e-03 | 0.0088456017 |
| g__Hungatella_A | 4.811802e-03 | 0.0095549574 |
| g__MGBC136627 | 4.364796e-03 | 0.0163003355 |
| g__Escherichia | 4.188365e-03 | 0.0266100578 |
| g__UMGS1251 | 4.159842e-03 | 0.0072716746 |
| g__Hungatella | 4.116388e-03 | 0.0190788279 |
| g__Clostridium_Q | 4.011996e-03 | 0.0052127439 |
| g__Brevinema | 4.009680e-03 | 0.0123080282 |
| g__Thomasclavelia | 3.902580e-03 | 0.0109042017 |
| g__Scatousia | 3.649941e-03 | 0.0102727788 |
| g__Enterocloster | 3.620384e-03 | 0.0047339282 |
| g__Mailhella | 3.612079e-03 | 0.0102470769 |
| g__Copromonas | 3.599368e-03 | 0.0050180473 |
| g__Ventrimonas | 3.513355e-03 | 0.0071233140 |
| g__Caccovivens | 3.343316e-03 | 0.0122194941 |
| g__Lawsonia | 3.289015e-03 | 0.0117180962 |
| g__Fournierella | 3.217226e-03 | 0.0062309921 |
| g__Limenecus | 3.163279e-03 | 0.0065731820 |
| g__MGBC133411 | 3.042418e-03 | 0.0074576633 |
| g__Mucinivorans | 2.900095e-03 | 0.0373193950 |
| g__Sarcina | 2.842205e-03 | 0.0171339212 |
| g__Acetatifactor | 2.738138e-03 | 0.0058449483 |
| g__Eisenbergiella | 2.697412e-03 | 0.0068331645 |
| g__Bacteroides_G | 2.682722e-03 | 0.0346150378 |
| g__CAJLXD01 | 2.633818e-03 | 0.0087495783 |
| g__Blautia | 2.558708e-03 | 0.0061407647 |
| g__C-19 | 2.275515e-03 | 0.0048691968 |
| g__Butyricimonas | 2.217402e-03 | 0.0050081070 |
| g__Velocimicrobium | 2.201224e-03 | 0.0066764022 |
| g__CAZU01 | 2.196385e-03 | 0.0065944622 |
| g__MGBC116941 | 2.160169e-03 | 0.0093827415 |
| g__Intestinimonas | 2.081898e-03 | 0.0039382613 |
| g__Negativibacillus | 2.069077e-03 | 0.0055137501 |
| g__Rikenella | 1.985389e-03 | 0.0037067264 |
| g__Phascolarctobacterium | 1.943804e-03 | 0.0050321931 |
| g__RGIG6463 | 1.789843e-03 | 0.0039682804 |
| g__JALFVM01 | 1.742360e-03 | 0.0038623852 |
| g__Oscillibacter | 1.510459e-03 | 0.0024992661 |
| g__Pseudoflavonifractor | 1.502976e-03 | 0.0027706264 |
| g__Acinetobacter | 1.485415e-03 | 0.0097446552 |
| g__Citrobacter_A | 1.392923e-03 | 0.0060348874 |
| g__Staphylococcus | 1.361709e-03 | 0.0050863070 |
| g__RGIG4733 | 1.303790e-03 | 0.0040480788 |
| g__UBA1436 | 1.214079e-03 | 0.0065017516 |
| g__Lachnotalea | 1.208197e-03 | 0.0049169585 |
| g__Ruthenibacterium | 1.202021e-03 | 0.0053632179 |
| g__14-2 | 1.184643e-03 | 0.0096299356 |
| g__Beduini | 1.173950e-03 | 0.0025142208 |
| g__Scatocola | 1.121193e-03 | 0.0044975564 |
| g__Faecisoma | 1.085585e-03 | 0.0055596648 |
| g__Enterococcus_A | 1.083991e-03 | 0.0099150522 |
| g__Faecimonas | 9.880526e-04 | 0.0083079244 |
| g__RGIG9287 | 9.638826e-04 | 0.0093228400 |
| g__CAG-345 | 9.490865e-04 | 0.0060198775 |
| g__Blautia_A | 9.208028e-04 | 0.0029082304 |
| g__RGIG1896 | 8.343610e-04 | 0.0062772562 |
| g__CAG-269 | 7.982448e-04 | 0.0047176841 |
| g__Marseille-P3106 | 7.938268e-04 | 0.0017618833 |
| g__WRHT01 | 6.429563e-04 | 0.0026979819 |
| g__Eggerthella | 6.407882e-04 | 0.0021266467 |
| g__IOR16 | 6.398942e-04 | 0.0020651222 |
| g__Anaerotruncus | 6.265445e-04 | 0.0016948517 |
| g__RUG14156 | 6.151780e-04 | 0.0022147136 |
| g__CHH4-2 | 6.145042e-04 | 0.0019997615 |
| g__Corynebacterium | 5.949530e-04 | 0.0060400054 |
| g__Serratia_A | 5.860616e-04 | 0.0076188009 |
| g__CAG-56 | 4.915432e-04 | 0.0016341973 |
| g__Merdimorpha | 4.736030e-04 | 0.0013057353 |
| g__MGBC140009 | 4.679334e-04 | 0.0024067320 |
| g__CALURL01 | 4.635287e-04 | 0.0016737486 |
| g__Aliarcobacter | 4.604457e-04 | 0.0050324221 |
| g__Scatenecus | 4.520215e-04 | 0.0019760876 |
| g__RGIG8482 | 4.399064e-04 | 0.0029753981 |
| g__Enterobacter | 4.073437e-04 | 0.0041317730 |
| g__Klebsiella | 4.054439e-04 | 0.0048910857 |
| g__Caccenecus | 3.941198e-04 | 0.0017802371 |
| g__Alcaligenes | 3.699430e-04 | 0.0048092594 |
| g__Plesiomonas | 3.633249e-04 | 0.0027105056 |
| g__HGM05232 | 3.617894e-04 | 0.0009776409 |
| g__JAHHSE01 | 3.592120e-04 | 0.0014900895 |
| g__Egerieousia | 3.419549e-04 | 0.0011583178 |
| g__Emergencia | 3.413630e-04 | 0.0017450341 |
| g__Enterococcus_B | 3.352316e-04 | 0.0022266910 |
| g__Stoquefichus | 3.026072e-04 | 0.0020503969 |
| g__Hepatoplasma | 2.989107e-04 | 0.0038858387 |
| g__Paracoccus | 2.959092e-04 | 0.0038468195 |
| g__Moheibacter | 2.771627e-04 | 0.0031476408 |
| g__Scatomorpha | 2.641015e-04 | 0.0010184339 |
| g__UBA7185 | 2.434328e-04 | 0.0014558191 |
| g__Eubacterium | 1.646822e-04 | 0.0006729067 |
| g__Sphingobacterium | 1.505775e-04 | 0.0012460017 |
| g__Devosia | 1.489995e-04 | 0.0015094318 |
| g__Anaerosporobacter | 1.454112e-04 | 0.0012747262 |
| g__Caccomorpha | 1.383122e-04 | 0.0010540604 |
| g__UBA2658 | 1.307451e-04 | 0.0007204965 |
| g__Protoclostridium | 1.277418e-04 | 0.0007646754 |
| g__Angelakisella | 1.268479e-04 | 0.0009221501 |
| g__Cetobacterium_A | 9.310217e-05 | 0.0008718575 |
| g__Rahnella | 6.470705e-05 | 0.0008411917 |
| g__Peptostreptococcus | 2.287338e-05 | 0.0002973539 |
genus_arrange <- genus_summary %>%
group_by(genus) %>%
summarise(mean=sum(relabun)) %>%
filter(genus != "g__")%>%
arrange(-mean) %>%
select(genus) %>%
mutate(genus= sub("^g__", "", genus)) %>%
pull()
genus_summary %>%
left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
mutate(genus= sub("^g__", "", genus)) %>%
filter(genus %in% genus_arrange[1:20]) %>%
mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")
4.3 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample)) %>%
full_join(functional, by = join_by(sample == sample))4.3.1 Wild samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric ) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.2 Acclimation samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.3 Antibiotics samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="2_Antibiotics") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.4 Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="3_Transplant1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.5 Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="4_Transplant2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.6 Post-Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.3.7 Post-Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)
4.4 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, tree = genome_tree)
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, dist = dist)4.5 Permanovas
4.5.1 1. Are the wild populations similar?
4.5.1.1 Wild: P.muralis vs P.liolepis
wild <- meta %>%
filter(time_point == "0_Wild")
# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))
wild_nmds <- sample_metadata %>%
filter(time_point == "0_Wild")4.5.1.3 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.631953 | 0.207198 | 6.795072 | 0.001 |
| Residual | 26 | 6.244347 | 0.792802 | NA | NA |
| Total | 27 | 7.876300 | 1.000000 | NA | NA |
4.5.1.4 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 2.018797 | 0.2566342 | 8.976049 | 0.001 |
| Residual | 26 | 5.847641 | 0.7433658 | NA | NA |
| Total | 27 | 7.866438 | 1.0000000 | NA | NA |
4.5.1.5 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.3786052 | 0.2108638 | 6.947419 | 0.001 |
| Residual | 26 | 1.4168908 | 0.7891362 | NA | NA |
| Total | 27 | 1.7954960 | 1.0000000 | NA | NA |
4.5.1.6 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0787916 | 0.1594096 | 4.930642 | 0.051 |
| Residual | 26 | 0.4154800 | 0.8405904 | NA | NA |
| Total | 27 | 0.4942716 | 1.0000000 | NA | NA |
Run 0 stress 0.1097679
Run 1 stress 0.1023647
... New best solution
... Procrustes: rmse 0.03374888 max resid 0.1642954
Run 2 stress 0.1244398
Run 3 stress 0.09968334
... New best solution
... Procrustes: rmse 0.01457553 max resid 0.06935865
Run 4 stress 0.09968336
... Procrustes: rmse 1.975136e-05 max resid 7.011925e-05
... Similar to previous best
Run 5 stress 0.09968328
... New best solution
... Procrustes: rmse 7.900375e-05 max resid 0.000299581
... Similar to previous best
Run 6 stress 0.1146655
Run 7 stress 0.09968333
... Procrustes: rmse 5.868544e-05 max resid 0.000219313
... Similar to previous best
Run 8 stress 0.09968376
... Procrustes: rmse 0.0002576756 max resid 0.0009812402
... Similar to previous best
Run 9 stress 0.1330762
Run 10 stress 0.1545632
Run 11 stress 0.1103823
Run 12 stress 0.09968331
... Procrustes: rmse 4.903019e-05 max resid 0.0001789655
... Similar to previous best
Run 13 stress 0.0996833
... Procrustes: rmse 4.110386e-05 max resid 0.0001546972
... Similar to previous best
Run 14 stress 0.09968344
... Procrustes: rmse 0.0001462737 max resid 0.0005606005
... Similar to previous best
Run 15 stress 0.1023647
Run 16 stress 0.1354406
Run 17 stress 0.1128821
Run 18 stress 0.09968348
... Procrustes: rmse 0.0002544292 max resid 0.0009536273
... Similar to previous best
Run 19 stress 0.1097682
Run 20 stress 0.1157848
*** Best solution repeated 7 timesRun 0 stress 0.009648577
Run 1 stress 0.01430506
Run 2 stress 0.01015941
Run 3 stress 0.01585762
Run 4 stress 0.01619152
Run 5 stress 0.01237722
Run 6 stress 0.01220831
Run 7 stress 0.0121089
Run 8 stress 0.01603707
Run 9 stress 0.009652863
... Procrustes: rmse 0.0238468 max resid 0.0879788
Run 10 stress 0.01454273
Run 11 stress 0.01641299
Run 12 stress 0.01428731
Run 13 stress 0.01201338
Run 14 stress 0.01202111
Run 15 stress 0.01613131
Run 16 stress 0.01433374
Run 17 stress 0.009835439
... Procrustes: rmse 0.005961207 max resid 0.02243503
Run 18 stress 0.01849256
Run 19 stress 0.009587194
... New best solution
... Procrustes: rmse 0.002722235 max resid 0.01055571
Run 20 stress 0.01191033
Run 21 stress 0.01231184
Run 22 stress 0.009481774
... New best solution
... Procrustes: rmse 0.01368353 max resid 0.05032182
Run 23 stress 0.01215389
Run 24 stress 0.01703314
Run 25 stress 0.01175194
Run 26 stress 0.01215852
Run 27 stress 0.01926281
Run 28 stress 0.01086729
Run 29 stress 0.01284959
Run 30 stress 0.009962763
... Procrustes: rmse 0.02558047 max resid 0.1009984
Run 31 stress 0.01267035
Run 32 stress 0.01897278
Run 33 stress 0.01028877
Run 34 stress 0.0103817
Run 35 stress 0.01431348
Run 36 stress 0.01202024
Run 37 stress 0.01497341
Run 38 stress 0.01211831
Run 39 stress 0.01172194
Run 40 stress 0.0116573
Run 41 stress 0.009543414
... Procrustes: rmse 0.01197212 max resid 0.04701637
Run 42 stress 0.01015838
Run 43 stress 0.01144335
Run 44 stress 0.01034566
Run 45 stress 0.01015069
Run 46 stress 0.009504455
... Procrustes: rmse 0.01042751 max resid 0.03821164
Run 47 stress 0.009896259
... Procrustes: rmse 0.0236055 max resid 0.09359012
Run 48 stress 0.01233124
Run 49 stress 0.01231938
Run 50 stress 0.01418098
Run 51 stress 0.01195001
Run 52 stress 0.009554307
... Procrustes: rmse 0.01304913 max resid 0.04807294
Run 53 stress 0.01209722
Run 54 stress 0.01737512
Run 55 stress 0.00972456
... Procrustes: rmse 0.0186762 max resid 0.07438517
Run 56 stress 0.01619304
Run 57 stress 0.009463878
... New best solution
... Procrustes: rmse 0.002238582 max resid 0.007837354
... Similar to previous best
*** Best solution repeated 1 times
4.5.2 2. Do the antibiotics work?
4.5.2.1 Acclimation vs antibiotics
treat <- meta %>% #antibiotics samples giving error when plotting
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))
treat_nmds <- sample_metadata %>%
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")4.5.2.2 Number of samples used
[1] 53beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)4.5.2.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.033529 | 0.0957379 | 6.818959 | 0.001 |
| species | 1 | 2.327685 | 0.1095866 | 7.805342 | 0.001 |
| individual | 26 | 9.722175 | 0.4577167 | 1.253885 | 0.003 |
| Residual | 24 | 7.157208 | 0.3369589 | NA | NA |
| Total | 52 | 21.240598 | 1.0000000 | NA | NA |
4.5.2.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.168854 | 0.1067496 | 9.686182 | 0.001 |
| species | 1 | 3.090152 | 0.1520953 | 13.800733 | 0.001 |
| individual | 26 | 9.684304 | 0.4766554 | 1.663479 | 0.001 |
| Residual | 24 | 5.373891 | 0.2644996 | NA | NA |
| Total | 52 | 20.317201 | 1.0000000 | NA | NA |
4.5.2.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.0671823 | 0.2208963 | 21.358337 | 0.001 |
| species | 1 | 0.8731487 | 0.0933035 | 9.021460 | 0.001 |
| individual | 26 | 4.0949687 | 0.4375828 | 1.627293 | 0.004 |
| Residual | 24 | 2.3228577 | 0.2482174 | NA | NA |
| Total | 52 | 9.3581574 | 1.0000000 | NA | NA |
4.5.2.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.1253769 | 0.4009905 | 40.2325710 | 0.001 |
| species | 1 | 0.0070451 | 0.0013292 | 0.1333606 | 0.777 |
| individual | 26 | 1.9000410 | 0.3584768 | 1.3833480 | 0.200 |
| Residual | 24 | 1.2678545 | 0.2392035 | NA | NA |
| Total | 52 | 5.3003175 | 1.0000000 | NA | NA |
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_treat <- beta_div_func_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
4.5.3 3. Do the FMT work?
4.5.3.1 Comparison between FMT2 vs Post-FMT2
transplant3 <- meta %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")
transplant3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts)),sort(as.character(rownames(transplant3))))
transplant3_nmds <- sample_metadata %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")4.5.3.2 Number of samples used
[1] 53beta_div_richness_transplant3<-hillpair(data=transplant3.counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3.counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3.counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3.counts, q=1, dist=dist)4.5.3.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.9809525 | 0.0778442 | 5.775712 | 0.001 |
| time_point | 1 | 0.7092107 | 0.0562799 | 4.175734 | 0.001 |
| type | 1 | 1.4479860 | 0.1149060 | 8.525540 | 0.001 |
| individual | 25 | 6.9157236 | 0.5488022 | 1.628753 | 0.001 |
| Residual | 15 | 2.5476146 | 0.2021678 | NA | NA |
| Total | 43 | 12.6014875 | 1.0000000 | NA | NA |
4.5.3.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.1101711 | 0.0895457 | 7.879642 | 0.001 |
| time_point | 1 | 0.7363935 | 0.0593971 | 5.226687 | 0.001 |
| type | 1 | 1.9027421 | 0.1534740 | 13.505059 | 0.001 |
| individual | 25 | 6.5351386 | 0.5271203 | 1.855374 | 0.001 |
| Residual | 15 | 2.1133659 | 0.1704628 | NA | NA |
| Total | 43 | 12.3978112 | 1.0000000 | NA | NA |
4.5.3.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.1341431 | 0.1005265 | 6.921553 | 0.001 |
| time_point | 1 | 0.1159136 | 0.0868654 | 5.980943 | 0.001 |
| type | 1 | 0.1423573 | 0.1066822 | 7.345390 | 0.001 |
| individual | 25 | 0.6512838 | 0.4880704 | 1.344205 | 0.070 |
| Residual | 15 | 0.2907075 | 0.2178554 | NA | NA |
| Total | 43 | 1.3344053 | 1.0000000 | NA | NA |
4.5.3.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | -0.0031096 | -0.0030953 | -0.1177247 | 0.813 |
| time_point | 1 | -0.0045864 | -0.0045653 | -0.1736359 | 0.840 |
| type | 1 | 0.0775554 | 0.0771987 | 2.9361453 | 0.141 |
| individual | 25 | 0.5385511 | 0.5360739 | 0.8155530 | 0.630 |
| Residual | 15 | 0.3962105 | 0.3943880 | NA | NA |
| Total | 43 | 1.0046211 | 1.0000000 | NA | NA |
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="top")
p1<-beta_neutral_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylo_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_func_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")
4.5.4 4. Are there differences between the control and the treatment group?
4.5.4.2 Number of samples used
[1] 27beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)4.5.4.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.6488906 | 0.0757032 | 2.115638 | 0.003 |
| type | 1 | 0.5615418 | 0.0655126 | 1.830847 | 0.010 |
| Residual | 24 | 7.3610760 | 0.8587842 | NA | NA |
| Total | 26 | 8.5715084 | 1.0000000 | NA | NA |
4.5.4.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.7984743 | 0.104299 | 3.065171 | 0.001 |
| type | 1 | 0.6051778 | 0.079050 | 2.323148 | 0.011 |
| Residual | 24 | 6.2519772 | 0.816651 | NA | NA |
| Total | 26 | 7.6556293 | 1.000000 | NA | NA |
4.5.4.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0700374 | 0.0635900 | 1.6594454 | 0.138 |
| type | 1 | 0.0184254 | 0.0167292 | 0.4365646 | 0.799 |
| Residual | 24 | 1.0129278 | 0.9196808 | NA | NA |
| Total | 26 | 1.1013906 | 1.0000000 | NA | NA |
4.5.4.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | -0.0048931 | -0.0058792 | -0.1628739 | 0.830 |
| type | 1 | 0.1161466 | 0.1395538 | 3.8660898 | 0.082 |
| Residual | 24 | 0.7210172 | 0.8663254 | NA | NA |
| Total | 26 | 0.8322707 | 1.0000000 | NA | NA |
p0<-beta_richness_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="top")
p1<-beta_neutral_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")
4.5.4.4 Number of samples used
[1] 28beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)4.5.4.4.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.8966107 | 0.1132233 | 3.515588 | 0.001 |
| type | 1 | 0.6463814 | 0.0816245 | 2.534445 | 0.001 |
| Residual | 25 | 6.3759661 | 0.8051521 | NA | NA |
| Total | 27 | 7.9189582 | 1.0000000 | NA | NA |
4.5.4.4.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.9398968 | 0.1224370 | 4.112304 | 0.001 |
| type | 1 | 1.0227481 | 0.1332297 | 4.474800 | 0.002 |
| Residual | 25 | 5.7139316 | 0.7443333 | NA | NA |
| Total | 27 | 7.6765766 | 1.0000000 | NA | NA |
4.5.4.4.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.1200820 | 0.1454923 | 4.647187 | 0.001 |
| type | 1 | 0.0592745 | 0.0718175 | 2.293931 | 0.029 |
| Residual | 25 | 0.6459931 | 0.7826902 | NA | NA |
| Total | 27 | 0.8253496 | 1.0000000 | NA | NA |
4.5.4.4.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0116738 | 0.0168395 | 0.4229082 | 0.526 |
| type | 1 | -0.0085273 | -0.0123007 | -0.3089208 | 0.917 |
| Residual | 25 | 0.6900904 | 0.9954612 | NA | NA |
| Total | 27 | 0.6932368 | 1.0000000 | NA | NA |
p0<-beta_richness_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="top")
p1<-beta_neutral_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")